/* Copyright 2011 Tobias Marschall, Email: tm@cwi.nl 
 * 
 * This file is part of string-cover.
 *
 * string-cover is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * string-cover is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with string-cover.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <iostream>
#include <limits>
#include <string.h>

#include "LagrangianBoundProvider.h"

using namespace std;

LagrangianBoundProvider::LagrangianBoundProvider(const Multipliers& initial_multipliers) : initial_multipliers(initial_multipliers) {
	// Step size parameter for Lagrangian optimization.
	this->initial_mu = 2.0;
	// Number of non-improving iterations after which the step size is decreased.
	this->allowed_nonimproving_iterations = 5;
	this->allowed_nonimproving_iterations_root = 30;
	// Factor by which mu is decreased.
	this->mu_decrease_factor = 0.5;
	// When mu drops below this threshold, we terminate.
	this->mu_threshold = 0.125;
	this->mu_threshold_root = 0.01;
	this->min_length = -1;
	this->max_length = -1;
	this->verbose = false;
}

LagrangianBoundProvider::~LagrangianBoundProvider() {
}

void LagrangianBoundProvider::setMinLength(int min_length) {
	this->min_length = min_length;
}

void LagrangianBoundProvider::setMaxLength(int max_length) {
	this->max_length = max_length;
}

void LagrangianBoundProvider::setVerbose(bool verbose) {
	this->verbose = verbose;
}

std::auto_ptr<Solution> LagrangianBoundProvider::compute_bounds(BranchAndBoundNode& node, double global_upper_bound, double* lower_bound, double* substring_weight_sums) const {
	// TODO: not nice: some code duplicated from SimpleBoundProvider
	const ProblemInstance& instance = node.getProblemInstance();
	const vector<string>& strings = instance.getStrings();
	const SubstringStore& substrings = instance.getSubstrings();
	const IntervalNumbering& intervals = instance.getIntervalNumbering();
#ifdef DEBUG
	cout << "Lagrange: global upper bound: "<< global_upper_bound << endl;
#endif
	Multipliers* multipliers;
	if (node.getMultipliers()==0) {
		multipliers = new Multipliers(initial_multipliers);
	} else {
		multipliers = new Multipliers(*node.getMultipliers());
	}
	// set all multipliers to zero for substrings included by this b&b node:
	// iterate over all intervals
	for (size_t string_nr=0; string_nr<strings.size(); ++string_nr) {
		for (size_t i=0; i<strings[string_nr].size(); ++i) {
			for (size_t j=i+1; j<=strings[string_nr].size(); ++j) {
				if (node.isIncluded(substrings.getIndex(string_nr, i, j))) {
					multipliers->set(string_nr, i, j, 0.0);
				}
			}
		}
	}
	double best_lower_bound = -numeric_limits<double>::infinity();
	// if all weights are integers, we can round up the lower bound
	double rounded_best_lower_bound = best_lower_bound;
	double current_lower_bound = -numeric_limits<double>::infinity();
	double best_upper_bound = numeric_limits<double>::infinity();
	auto_ptr<vector<StringFactorization> > best_factorizations;
	// we call an iteration "successful", if it make the gap smaller by at least 1%
	int nonsuccessful_iterations = 0;
	int nonsuccessful_iterations_with_same_mu = 0;
	// size of gap that resulted from last successful iteration
	double best_successful_gap = numeric_limits<double>::infinity();
	Multipliers* last_multipliers = 0;
	double allowed_iterations = node.isRoot()?allowed_nonimproving_iterations_root:allowed_nonimproving_iterations;
	// decrease mu until it drops under a threshold
	for (double mu=initial_mu; mu>=(node.isRoot()?mu_threshold_root:mu_threshold); mu*=mu_decrease_factor) {
		if (verbose) {
			cout << "Lagrange: mu="<< mu << endl;
		}
		nonsuccessful_iterations_with_same_mu = 0;
		// iterate until a certain number of non-improving iterations have been done
		// for the current mu
		int nonimproving_iterations = 0;
		while (true) {
#ifdef DEBUG
			cout << "Lagrange: using these multipliers:"<< endl << *multipliers;
#endif
			// bitset of selected x-variables
			boost::dynamic_bitset<> selected_substrings(substrings.size());
			// bitset of selected z-variables
			boost::dynamic_bitset<> selected_intervals(intervals.size());
			if (substring_weight_sums != 0) {
				memset(substring_weight_sums, 0, sizeof(double)*substrings.size());
			}
			auto_ptr<vector<StringFactorization> > factorizations(new vector<StringFactorization>());
			current_lower_bound = 0.0;
			// determine which x-variables need to be set to one:
			//   1) compute the coefficients: w(t) - \sum_{s[i..j]=t}\lambda_{s,i,j}
			double* x_coefficients = new double[substrings.size()];
			for (size_t i=0; i<substrings.size(); ++i) {
				x_coefficients[i] = instance.getWeight(i);
			}
			// iterate over all intervals
			for (size_t string_nr=0; string_nr<strings.size(); ++string_nr) {
				for (size_t i=0; i<strings[string_nr].size(); ++i) {
					for (size_t j=i+1; j<=strings[string_nr].size(); ++j) {
						x_coefficients[substrings.getIndex(string_nr, i, j)] -= multipliers->get(string_nr, i, j);
					}
				}
			}
			//   2) use all x-variables that have a negative coefficient
			for (size_t i=0; i<substrings.size(); ++i) {
				if (node.isForbidden(i))	continue;
				if (node.isIncluded(i) || (x_coefficients[i]<0.0)) {
					selected_substrings.set(i);
					current_lower_bound += x_coefficients[i];
				}
			}
			delete[] x_coefficients;
			// compute shortest path
			boost::dynamic_bitset<> feasible_solution(substrings.size());
			for (size_t string_idx=0; string_idx<strings.size(); ++string_idx) {
				const string& s = strings[string_idx];
				// min_costs[i] is the shortest path length to the point between s[i-1] and s[i]
				double* min_costs = new double[s.size()+1];
				int* back_trace = new int[s.size()+1];
				min_costs[0] = 0.0;
				back_trace[0] = -1;
				for (size_t i=1; i<=s.size(); ++i) {
					min_costs[i] = numeric_limits<double>::infinity();
					back_trace[i] = -1;
					for (size_t j=0; j<i; ++j) {
						// index of substring for edge j->i
						size_t substring_idx = substrings.getIndex(string_idx,j,i);
						int length = i-j;
						assert(length>0);
						if (node.isForbidden(substring_idx)) continue;
						if ((min_length>0) && (length<min_length)) continue;
						if ((max_length>0) && (length>max_length)) continue;
						double edge_weight;
						if (node.isIncluded(substring_idx)) {
							edge_weight = 0.0;
						} else {
							edge_weight = multipliers->get(string_idx, j, i);
						}
						double c = min_costs[j] + edge_weight;
						if (c<min_costs[i]) {
							min_costs[i] = c;
							back_trace[i] = j;
						}
					}
				}
				if (min_costs[s.size()] < numeric_limits<double>::infinity()) {
					current_lower_bound += min_costs[s.size()];
				} else {
#ifdef DEBUG
					cout << "Lagrange: no feasible solution exists in this b&b node." << endl;
#endif
					delete[] back_trace;
					delete[] min_costs;
					delete multipliers;
					return auto_ptr<Solution>(0);
				}
				// use backtrace to obtain factorization
				StringFactorization factorization(s.size());
				int i = s.size();
				while (back_trace[i]>=0) {
					size_t substring_nr = substrings.getIndex(string_idx,back_trace[i],i);
					if (substring_weight_sums != 0) {
						substring_weight_sums[substring_nr] += multipliers->get(string_idx,back_trace[i],i);
					}
					feasible_solution.set(substring_nr);
					selected_intervals.set(intervals.getIndex(string_idx,back_trace[i],i));
					if (back_trace[i]==0) break;
					i = back_trace[i];
					factorization.addBreakpoint(i-1);
				}
				factorizations->push_back(factorization);
				delete[] back_trace;
				delete[] min_costs;
			}
			if ((last_multipliers!=0) && (last_multipliers!=multipliers)) delete last_multipliers;
			last_multipliers = multipliers;
			// compute upper_bound: score of current feasible solution
			double upper_bound = 0.0;
			size_t n = feasible_solution.find_first();
			for (; n!=boost::dynamic_bitset<>::npos; n=feasible_solution.find_next(n)) {
				upper_bound += instance.getWeight(n);
			}
#ifdef DEBUG
			cout << "Lagrange: factorizations for feasible solution (score "<< upper_bound <<"):" << endl;
			for (size_t i=0; i<strings.size(); ++i) {
				cout << factorizations->at(i).printFactorizedString(strings[i]) << endl;
			}
#endif
			if (current_lower_bound>best_lower_bound) {
#ifdef DEBUG
				cout << "Lagrange: Improved lower bound from " << best_lower_bound << " to " << current_lower_bound << endl;
#endif
				best_lower_bound = current_lower_bound;
				if (node.getInstance().hasIntegerWeights()) {
					rounded_best_lower_bound = ceil(best_lower_bound-1e-10);
				} else {
					rounded_best_lower_bound = best_lower_bound;
				}
				nonimproving_iterations = 0;
			} else {
				nonimproving_iterations += 1;
			}
			if (upper_bound < best_upper_bound) {
				best_factorizations = factorizations;
#ifdef DEBUG
				cout << "Lagrange: Improved upper bound from " << best_upper_bound << " to " << upper_bound << endl;
#endif
				best_upper_bound = upper_bound;
				nonimproving_iterations = 0;
			}
			double current_gap = 1.0-current_lower_bound/best_upper_bound;
			if (current_gap/best_successful_gap <= 0.99) {
#ifdef DEBUG
				cout << "Lagrange: Improved best successful gap from " << best_successful_gap << " to " << current_gap << endl;
#endif
				nonsuccessful_iterations = 0;
				nonsuccessful_iterations_with_same_mu = 0;
				best_successful_gap = current_gap;
			} else {
				nonsuccessful_iterations += 1;
				nonsuccessful_iterations_with_same_mu += 1;
			}
			if (verbose) {
				cout << "Lagrage: resulting lower bound:"<< current_lower_bound << ", 1.0-lower/upper = 1.0-" << current_lower_bound << '/' << best_upper_bound << " = " << current_gap << " (" << nonimproving_iterations << '/' << nonsuccessful_iterations_with_same_mu << '/' << nonsuccessful_iterations << ")"<< endl;
			}
			if (rounded_best_lower_bound>=global_upper_bound) break;
			if (nonimproving_iterations > allowed_iterations) break;
			if (nonsuccessful_iterations_with_same_mu > 3*allowed_iterations) break;
			if (nonsuccessful_iterations > 5*allowed_iterations) break;
			// quit if optimum has been reached
			if (fabs(1.0-rounded_best_lower_bound/best_upper_bound) < 1e-10) break;
			// otherwise adapt multipliers
			multipliers = multipliers->adapt(mu, min(best_upper_bound, global_upper_bound), current_lower_bound, selected_substrings, selected_intervals).release();
			// quit if all gradients have been zero
			if (multipliers==0) break;
		}
		if (rounded_best_lower_bound>=global_upper_bound) break;
		if (nonsuccessful_iterations > 5*allowed_iterations) break;
		if (fabs(1.0-rounded_best_lower_bound/best_upper_bound) < 1e-10) break;
		if (multipliers==0) break;
	}
	if (substring_weight_sums != 0) {
		double* multiplier_sums = new double[substrings.size()];
		last_multipliers->multiplierSums(multiplier_sums);
		for (size_t i=0; i<substrings.size(); ++i) {
			substring_weight_sums[i] /= multiplier_sums[i];
		}
		delete[] multiplier_sums;
	}
	*lower_bound = rounded_best_lower_bound;
	if ((last_multipliers!=0) && (last_multipliers!=multipliers)) delete last_multipliers;
	node.setMultipliers(auto_ptr<Multipliers>(multipliers));
	return auto_ptr<Solution>(new Solution(instance,best_factorizations));
}

void LagrangianBoundProvider::setParameters(double initial_mu, int allowed_nonimproving_iterations, int allowed_nonimproving_iterations_root, double mu_decrease_factor, double mu_threshold, double mu_threshold_root) {
	this->initial_mu = initial_mu;
	this->allowed_nonimproving_iterations = allowed_nonimproving_iterations;
	this->allowed_nonimproving_iterations_root = allowed_nonimproving_iterations_root;
	this->mu_decrease_factor = mu_decrease_factor;
	this->mu_threshold = mu_threshold;
	this->mu_threshold_root = mu_threshold_root;
}


